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PREFACE 


The purpose of this final report prepared under NASA Contract 
NAS 5-21680 is to summarize the work carried out during the term of 
this contract. Two sets of programs were developed for computing 
characteristics of the scattered radiation emerging at selected levels 
of a plane-parallel, nonhomogeneous atmosphere containing an 
arbitrary vertical distribution of ozone concentration and/or aerosol 
(spherical particles of known refractive index) number-density, and 
bounded at the lower end by a Lambert ground of known reflectivity. 

The atmosphere is assumed to be illuminated at the top by unidirectional, 
unpolarized, monochromatic radiation. 

In the first set of programs referred to as a scalar case, polarization 
characteristics of the scattered radiation are neglected during evaluation 
of the contribution to the emergent radiation from higher orders of 
scattering. Consequently, this first set of programs enables the user 
to compute, with moderate amount of computer time and storage require- 
ments, variations of the scattered intensity as a function of several 
parameters such as, zenith angle, azimuth angle, Lambert ground re- 
flectivity, sun's position, wavelength, and atmospheric parameters. 

The second set (referred to as a vector case) of programs can be used 
to compute all four characteristics (intensity, as well as degree, 
direction and ellipticity of polarization) of the scattered radiation for the 
various parameters listed above. 

Each set consists of four separate programs. The first two 
programs are used for evaluating needed scattering properties of a 
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unit volume containing a known size distribution of spherical particles 
made of a material whose refractive index with respect to air is given. 
The remaining two programs are for evaluating various characteristics 
of the scattered radiation by using a modified Fourier series method 
combined with the Gauss-Seidel Iterative procedure. 
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I. 


INTRODUCTION 


The analysis of data acquired with the Backscatter Ultraviolet 
(BUV) experiment aboard the Nimbus IV satellite requires knowledge 
of the respective contributions to the outgoing radiation due to 
scattering by molecules and aerosols, as well as due to reflection 
by the surface underlying the atmosphere. The objective of NASA 
Contract No. NAS 5-21680 is to develop and to test programs aimed 
at computing various characteristics of the radiation emerging at 
selected levels of a plane-parallel, nonhomogeneous atmosphere 
containing an arbitrary vertical distribution of ozone and/or aerosol 
number-density, and bounded at the lower end by a Lambert ground 
of known reflectivity. The scattering aerosols are assumed to be 
spherical particles made up of a material whose refractive index with 
respect to air is known. It is further assumed that this refractive 
index and the size-distribut.'on function are independent of height. 

The atmosphere is assumed ! .o be illuminated at the top by unidirectional, 
unpolarizad , monochromatic radiation. Various characteristics of the 
scattered radiation emerging at selected levels of a given atmospheric 
model are evaluated by using a modified Fourier series method {Dave 
and Gazdag, 1970, Dave, 1970A and 1970B) combined with the Gauss- 
Seidel iterative procedure (Hildebrand, 1956, Chapter 10). 

The computer requirements for evaluating scattered-radiation 
characteristics depend upon several factors such as optical properties 
of the model, wavelength to be investigated, desired accuracy, and the 
number of cases for which output is required for a meaningful analysis 
of a given problem. The total computer requirement for a given study 


can heavily tax some of the most modern computing facilities available 
to typical research group associated with a large organization. Hence, 
we have developed, and tested, two sets of programs which are described 
below. 

In the first set of programs referred to as a scalar case, polarization 
characteristics of the scattered radiation are neglected during evaluation 
of the contribution to the emergent radiation from the radiation which Is 
scattered more than once in the atmosphere. Consequently, this first 
set of programs enables the user to evaluate, with moderate amounts of 
computer time and storage requirements, variations of the scattered 
intensity as a function of several parameters such as zenith (or nadir) 
angle, azimuth angle, Lambert ground reflectivity, sun's position, 
wavelength, and atmospheric parameters. It is known that the values of 
intensities obtained with such a scalar treatment of the transfer problem 
can differ significantly {~ 10 to 30%) from those obtained with a full 
vector treatment, i.e. , when polarization of the scattered radiation is 
properly accounted for (e.g., Chandrasekhar, 1950, Sec. 70.4; Adams 
and Kattawar, 1970; Dave and Gazdag, 1970). The magnitude of the 
difference depends upon the optical thickness of the atmosphere, 
zenith angle of the sun and direction of observation, location of the 
observer within the atmosphere, and nature of the scattering phase 
function. However, because of the excessive computer requirements , 
some investigators not concerned with polarization aspects of the 
scattered radiation work with numerical results obtained with such a 
scalar program. A plausible reason advanced for justifying the validity 
of the ultimate findings of such an investigation is that of biasing of 
all numerical results by the same amount. The second set of programs 
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can be used to compute all four characteristics (intensity, as well as 
degree, direction and ellipticlty of polarization) of the scattered 
radiation for the various , irameters listed above. 


Each set of the programs described above consists of four 
different programs. A description of the purpose, computational 
procedure, FORTRAN listing, and the results of test runs obtained 
with each of these eight programs are given in the following eight 
reports prepared for this contract: 

Title of the Report Date Submitted No, of Pages 

Technical Report - Scalar Case, Program I 20 Ian. 1972 43 

Technical Report - Scalar Case, Program II 24 Feb. 1972 38 

Technical Report - Scalar Case, Program III 06 Mar. 1972 50 

Technical Report - Scalar Case, Program IV 04 May 1972 106 

Technical Report - Vector Case, Program I 10 Feb. 1972 50 

Technical Report - Vector Case, Program II 18 May 1972 46 

Technical Report - Vector Case, Program III 08 June 1972 70 

Technical Report - Vector Case, Program IV 03 August 1972 138 

A brief description of the basic radiative transfer equation used 
for this work, and the function of each of the eight programs listed 
above are given in the sections which follow. 


II. DEFINITION Of THE BASIC QUANTITIES 



In this section we shall give brief definitions of various quantities 
required for writing down the eq lation of radiative transfer for a plane- 
parallel, nonhomogeneous atmosphere. 


2.1 Directional parameters The direction of propagation 

is denoted by the zenith angle a ( = cos Mul ) which the direction 
makes with local zenith, and the azimuth angle o which the 
meridian plane containing the direction of propagation makes with 
an arbltrari.y chosen meridian plane. A negative (positive) sign 
before n will imply the radiation traveling along the downward 
(upward) direction. A subscript “o'* for and c p will imply the 
direction of propagation of the incident solar radiation. 


2.2 Normal Rayleigh- scattering optical depth t ' : This 
optical depth of a level situated at a height of hkm above the 
ground is given by 


T (s ' r) = | k (s) p(h) dh , (1) 

h 

where 

( s) 

k = mass scattering coefficient of air at wavelength X , 

and 

p (h) = density of air at height h . 
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(s m) 

2.3 Normal Mle-scattertng optical depth t * : The attenuation 

suffered by a beam of radiation travelling along the vertical due 
to scattering by atmospheric aerosols is determined after making 
the following assumptions: (a) aerosol particles are spherical in 
shape, (b) all particles are made of «.ne material whose refractive 
index with respect to air is mf^nj-in^, and (c) the size-distri- 
bution function, l.e. , n(r) ® number of particles per cc in 1 micron 
interval at adius r , is independent of height. Accordingly, we 
have , ! 

r (s ' m) = 3 N(h)/N (2) 

s o 

where 

0 = volume scattering coefficient per cm [Eq. (3)], 
s 

N(h) - total number of particles in one sq. cm column above 

the level located at a height of h km above the ground, 
and 

N q - total number of particles per cc [Eq. (4) ] . 

The explicit expression for 0 is 

3 

r 

f max 2 

0 = u l 0 (x,m) r n(r) dr , (3) 

s j s 

r min 

where 

Q (x,m) = efficiency factor for scattering for a particle with 
s 

size parameter x = 2nr/X ; r is the radius of the 
particle (Van de Hulst, 1957, Sec. 9.3,2). 
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The quantity N q is given by 
r 

- max 

N o " \ dr ' (4) 

r min 

f a m) 

2.4 Normal Mle absorption op tica l depth r Va, ‘‘ : Following 
the procedure given in Sec. 2.3, the Mie absorption optical 
depth of a level is given by 

r (a,m) = a N(h)/N (5) 

a o 

where the explicit expression for 6 is obtained after replacing 

a 

s with a in Eq. (4); Q (x,m) is the efficiency factor for 

a 

absorption. 

(a) 

2.5 Normal ozone-absorption optical depth r : This quantity 
is given by 

= a x(h) , (6) 

where 

a = absorption coefficient (to the base e , e.g. , Handbook 
of Geophysics, Revised Edition, U.S. Air Force, 
Macmillan Company, 1961, pp. 16-23 to 16-25) of 
one cm-atm of ozone at NTP, for the wavelength X , and 

x(h) ~ total amount of ozone in cm-atm above the level. 
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2.6 Normal total optical depth r : The attenuation suffered 

by a beam of radiation travelling (in a direction making an angle 
8 with local zenith) from the top of the atmosphere to a level 
located at a height h km above the ground is exp(-r/V) • The 
quantity r is given by 


t ~ T 


(r) , _(m) 


+ r 


( 7 ) 


where the attenuation optical depth due to gases is given by 


.(r) _ (s,r) (a) 

— T -r T 


( 8 ) 


and the attenuation optical depth due to aerosols is given by 


(m) (s,m) . (a,m) 

T — T *«- T 


(9) 


2,7 Albedo of single scattering ^(t) : This quantity defines 

the relationship between absorption and scattering characteristics 
of the materials located at the level r . It is given by 


oj(t) - 


a T (s - m) + & T (s ' r) 

. (m) (r) 

Ar + At 


( 10 ) 


where At^^ represents a change in the -th optical depth. 


2.8 Turbidity factor T(r) : This quantity defines the relationship 
between the scattering characteristics of molecules and aerosols 
located at the level r . It is defined as follows: 
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T(t) = 


( 11 ) 


2.9 Stokes parameters J ^r;u,o) : Since we are Interested in 
all aspects of scattering, it is necessary to treat a beam of 
monochromatic radiation of wavelength > as a vector or as a 
one-column matrix with four elements. In the Stokes representation 
(Chandrasekhar, 1950, Sec. 15), the intensity matrix is represented 
by 

Kt;m,0) = 


I g (t;m,0) 

I u (T;*i.«P) 
I v (T;n,cp) 


( 12 ) 


The first two elements [I (T;^,^) and I(T;ji,<p)] of this matrix 

e r 

represent the specific intensities (or for short, intensities) of the 
beam in two directions , e and r , respectively. These two 
directions are mutually at right angles to each other, such that 
the e-r plane is perpendicular to the direction of propagation of 
energy. The intensity of the beam, i.e. , the amount of energy 
per unit wavelength interval at k flowing per unit time in a cone 
of unit solid angle with axis in the direction n,<p , is given by 


I(t;h,0) = I + I (t;m,<p) . (13) 

e r 

The degree of polarization of the beam is given by 
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1/2 

P - [0 2 (t;m.<p) + 1 2 (t;u,<p) + I 2 (r;u,tp)] /I(t :n,tp) , (14) 


where 


Q{t;h,v) = I 3 (t;m,<p) - I (T;n,<p) . 

V I 


(15) 


The angle (X) which the direction of maximum vibration 
makes with the e-direction, i.e. , the deviation of the plane of 
polarization, or alternately the direction of polarization, is 
given by 


tan 2 x = I u (r;fx,<p)/ • 


(16) 


The ellipticity of polarization defined as the ratio of minor 
to the major axis of the ellipse traced by the electric vector, is 
given by 


2 2 2 

tan 8 = - I v (r;p,<p)/ [Q (t;m,<p) + I u (T;*i,<p) + I v (r;^,«p)] 


1/2 


2 2 
+ [Q (r;y,<p) + I u (r;M,<p)] 


1/2 


(17) 


It may be noted that the quantities P, X , and 6 given by 
Eqs. (14), (16), and (17) as well as the intensity I(T;fX,<p) 
given by Eq. (13) are functions of directional parameters 

(-^ 0 , <P Q ) of the incident radiation, and various optical properties 
of the atmosphere and underlying surface. 
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2.10 Normalized scattering phase matrix P(t;u , o; u' ,o') : 

This four-by-four matrix gives, on a relative scale, the Stokes 
parameters of the radiation scattered by a unit volume at the 
level t in the direction when the volume is illuminated 

by an incident beam of known Stokes characteristics from the 
direction m' , «P' * For a volume containing both aerosols and 
molecules, this matrix has the following form: 

P n' ,<p‘) = T (r) M(m,<p; n' , 0 ') + [l-T(r)] R(m,o:m‘ 

(18) 

where T(r) is the turbidity factor defined above. The quantities 
M(m,<P;m' ,<P') and M' . <P ') are four-by-four normalized 

matrixes for a unit volume of aerosol (Mie scattering) and air 
(Rayleigh scattering), respectively. Their ij-th elements are 
either even or odd functions of ip' - <p depending upon the value of 
the subscript ij (Dave, 1970A). The azimuth dependence is 
given by: 

M^(fi,«p; n' ,<p') = y] m|^(u,m') * (19) 

n = 1 

and 

3 

R^(ji,<p; n' ,<p') « ~ ^ 1 ( 2 °) 

n= 1 

where for ij = 11, 12, 21, 22, 33, 34, 43, and 44, these 
elements are even functions, i.e. , 

f . ( <p' - <p) = cos[(n-l) (<p' - <p)] , (21) 

n - 1 
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and for ij-13, 14, 23 , 2*5 , 31, 32 , 41, and 42, those elements 
are odd functions, i.e. , 

f (< 0'-ip)= sin ( (n— 1 ) ( tp ' **«P) 1 • ( 22 ) 

n- L 

(n) (n) 

The explicit expressions for m| (ji,y‘) and Rjj (M#M’) can be 
found in Technical Reports - Vector Case, Programs III and IV 
respectively. With the Fourier series form of the Mle and 
Rayleigh phase matrices given above, the ij ~th clement of the 
matrix P (r;^,o; n' ,o‘) can be written as 

N(jiiM') 

P^(r;M,tp ; n' ,(p') = pj^fri/i.^') f n _ 1 (o‘ -<p) , (23) 

n = 1 

where 

p[j l) (r;4 ,M‘) = T(t) + [ 1 - T(t) ] R^V^’l . (24) 

Evidently, r|' j (m,M') £ 0 for n>3 . The series given by 
Eq. (23) is even or odd depending upon the value of subscript ij . 
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III. EQUATION OF RADIATIVE TRANSFER 


The basic transfer equation for a plane-parallel atmosphere of 
infinite extent and homogeneous along the horizontal, but of finite 
extent and nonhomogeneous in the vertical has the following form. 

- ji(t) , (25) 

dT ■ r w 

where the source matrix J[t;h,<p) is given by 
1 - T/ "o 

I (t;m.p) = ;e p(r :n,<p; -u , <p ) • F 

H *-» 0 0 *"* 

+1 2rr 

+ 4^ f ^ P (r;n,(p; , 0 ') * Urif*' ,<p’) dp' d(fi’ . (20) 

-1 0 

The boundary conditions are given by 
,1(0; -n, (fl) = 0 , 

and , (27) 

I(t. 5 0 

D 

where r, is the total normal optical thickness of the atmosphere. In 
b 

arriving at Eq. (26) for the source matrix, it is assumed that the atmos- 
phere Is illuminated by unidirectional ( - M 0 < <P Q ) solar radiation with 
a net flux of ttF per unit area normal to the direction of propagation. 

The Stokes parameters of the incident solar radiation are represented 
by a one-column matrix given by 
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« r - 


TT 


{ 28 ) 


V 

o 
F 

r 

r 

u 
F 

v 

The net flux np - n(F + r ) . For our calculations wo have assumed 

e r 

that F = F = 1/2, and F - F = 0, i.e., the Incident solar Is 
e r u v 

unpolarized, and has a net flux *t per unit area normal to its direction 
of propagation. All the quantities appearing in Eqs. (25) to (27) are 
defined in Sec. 2. It may be noted that this transfer equation is, 
strictly speaking, valid only ior a monochromatic radiation. However, 
one may us: the same program for computing various characteristics 
of the emergent radiation integrated over a specified wavelength interval 
provided he feels justified in using average values of various input 
parameters described in Secs. 2.2 to 2.S. 

The form of the normalized phase matrix for a unit volume [Eq. (23)] 
suggests that we expand the elements of I(t;u,«p) and 1(t;m,<p) in 
similar forms. After substituting the Fourier series expressions for 
l(r;^,tp) , J(t;m-<p) and P(t;M/<P; ,<P') in Eqs. (25) and (26) , we 

can carry out the integration over azimuth angle (<p') analytically by 
making use of the well-known orthogonality relations obeyed by the 
trigonometric functions. It can be shown that (Technical Report - 
Vector Case, Program IV, Sec. 2.2) Eq. (25) reduces to the following 
system of uncoupled integro-differential equations: 
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(29) 


d . ,W ( 

dT 


(r) J^(r:^) 


where the superscript n increases from 1 to N(m q ) • The magnitude 
of this upper limit depends upon the nature of aerosols in the atmospheric 
model under investigation, and also on the zenith angle of the sun 
(Dave and Gazdag, 1970). The n-th term of the Fourier series for the 
source matrix has the following form: 



(Tfji) 


l 

4 



(t;m. -M 0 > 





+ l 


P (n) (r;MiM') • I {n) (r;M') d*i‘ 


(30) 


-1 


where 6, is the Kronecker delta function given by i. =1 when 
In In 

n = 1 and otherwise zero. The boundary conditions for the n-th term 
of Eq. (29) are given by 


and 


l (n) {0;-*i) s 0 


I (n) (T.;+M) « 0 

D 


(31) 


Identical equations for the scalar case are arrived at by treating 
the quantities I (r;/i,<p) , J(t;m.<p) * P(t;m,<p; h' ,<£>’) , I^(r;ji) * 

N*/" 

J^(t;k) . J^(r;u) , and P^(t;m»<p) as scalar quantities (Technical 
Report - Scalar Case, Program IV, Sec. 2.2). 
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In order to evaluate the contribution to 1 (t;u,c>) due to the 
presence of an idealised Lambert ground located underneath the atmos- 
phere, we can write equations similar to Cqs. (25) to (27) but with 
an isotropic source illuminating the model from below. In this case, 
we can easily carry out reductions similar to the one used for arriving 
at Eqs. (29) through (31), When this is done (Technical Report - Vector 
Case, Program IV, Sec. 2,3, or Technical Report - Scalar Case, Program 
IV, Sec. 2,3), we find that because of the nature of the illuminating 
source, this contribution to I(r;ji<o) is azimuth independent. 
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IV. FUNCTION or VARIOUS PROGRAMS 


The job of evaluating characteristics of the scattered radiation 

emerging at selected levels of a nonhomogeneous , plane-parallel 

atmosphere can be divided into the following two main tasks: (a) 

(n) 

Computations of the elements of the phase matrix P (T/p.p 1 ) given 
by Eq. (24), and (b) numerical solution of the transfer equation given 
by Eqs. (29) through (31) for the required number of terms of the Fourier 
series. As mentioned in Sec. I of this report, we have provided two 
sets of programs; one set of four programs for the scalar case, and 
another set again consisting of four programs for the vector case. 

In what follows, we will describe the functions of four programs for 
the vector case. The i-th program for the scalar case performs functions 
identical to the one performed by the i-th program for the vector case, 
but by treating related vector quantities as scalars. 

The function of the first program is to enable the user to compute 
values of coefficients of the Legendre series for four di.ierent scattering 
functions which can be used for evaluating the intensity as well as the 
degree, direction, and ellipticity of polarization of the radiation 
scattered by a sphere of known refractive index. Values of efficiency 
factors for extinction, scattering, and absorption for the sphere under 
investigation are also computed. This program enables the user to 
compute, and to store on a magnetic tape for further use with the second 
program, values of these quantities for hundreds of discrete but equally 
spaced values of the size parameter of the sphere. 
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The function of the second program is to enable the us*, i to compute 
coefficients of the Legendre series representing four scattering functions 
of the normalized scattering phase matrix of a unit volume containing a 
known size distribution of spherical particles made of the same material. 
Values of the volume extinction and volume scattering coefficients are 
also obtained during execution of this program. Values of these 
quantities are stored In punched-card form for use with the third program. 


The function of the third program is to enable the user to compute 

(n) 

the Fourier coefficients [appearing in Eqs. (19) and (24) ] 

of series representing elements of a normalized lour-by-four scattering 
phase matrix of a unit volume containing an arbitrary distribution of 
spherical particles. These coefficients are computed for the values of 
9' (-cos ^|fi' |) and 9 (-cos ^|fi() given by m 1 = 0°(2°) 180° , and 
9 = 0 J (2°) 90° . This program also prepares a table of the upper limit 
N(fi' ,*i) of the Fourier series [Eq. (19) ]. For further use with the 
fourth program , values of M|” (ji,,P) for all values of n, i, j, 9’ , and 
9 are stored on a magnetic tape, and those of N(m' ./il are stored in 
punched-card form. 


The function of the fourth program is to enable the user to compute 
the intensity (I) , the degree of polarization (P), the direction of 
polarization ( x ) » and the ellipticity of polarization (tan 3) of the 
scattered radiation emerging at selected levels of the atmospheric 
model under investigation. After completing computations for all re- 
quired Fourier components of the intensity [Eq. (29)1 values of I , P , 

X , and tan 3 of the scattered radiation emerging at selected levels 
of atmospheric models are computed and printed out for the desired 
values of Lambert ground reflectivities, and azimuth as well as zenith 
angles of the directions of the scattered radiation. 


17 


V. 


RELIABILITY 01' NUMERICAL RESULTS 


Because of the nature of the computations and also because of 
the unavailability of accurate numerical results in published form, it 
is not possible to make a definite statement about the absolute accuracy 
ot the numerical values obtained using these programs. However, based 
on various considerations and direct as well as indirect comparison with 
published and unpublished results of higher quality, we make the 
following statements: 

Various Legendre coefficients as well as values of the efficiency 
factors for absorption and scattering obtained using the first program for 
the vector {or scalar) case, are e V (.-’Cted to be correct to the first ten 
significant figures in most cases. 

The reliability of the normalized Legendre coefficients and of the 
scattering and absorption volume coefficients of a unit volume con- 
taining a known size distribution (analytic function) of spherical 
particles (second program for the Vector or the Scalar Case), is 
primarily determined by the value of the integration increment (Ax) 
used for the computations. For typical atmospheric aerosol distri- 
butions, the choice of Ax is governed to a considerable extent by 
the economics of computations. For such cases, a value of about 
0,1 for Ax is generally considered fine enough to provide results 
of about four significant figure reliability. 

Values of the Fourier coefficients of a series representing 
elements of a normalized four-by-four scattering phase matrix (or of 
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a series representing phase function for the scalar case) of a unit 
volume as obtained using the third program are expected to be as 
accurate as values of the normalized Legendre coefficients used as 
input. 


The reliability of the values of the intensity tor the scalar case 
and of all four parameters for the vector case is primarily determined 
by the value of the Integration increment ( A9 ) used for evaluating 
the contribution to the source term by atmospheric radiation [e.g. , the 
second term in Eq. (30)]. A provision is made in the fourth program 
to enable the user to work with a value of 2°, 6°, or 10° for A 8 . 

The reliability of this output also depends to some extent on the 
nature of the scattering phase matrix used in the model. Thus, for a 
given A 6 , the numerical results for an aerosol model with a very 
strong forward peak (~ 300 or more terms in Legendre representation) 
are bound to be less accurate than those for a model with a less strong 
forward peak. 

Ocher input parameters which determine the accuracy of the 
numerical output of the fourth program are the number of layers into 
which the basic atmospheric model is Ivided, the nature of vertical 
nonhomogeneities, and the value of the albedo of single scattering. 

A provision is made in the fourth program for permitting the user to 
divide his basic atmospheric model into 160/n number of layers where 
the integer n can assume anyone of the following values: 1, 2, 4, 

5, 8, or 10. The results obtained with finer division of the atmosphere 
can be expected to be more accurate than those obtained with the 
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coarser ones. However, the value ot the parameter ri and that of 
described In the preceding paragraph should be compatible. It is 
recommended that for ;. i > - 2 ° , the atmosphere be divided Into a 
sufficient number of layers so that the total scattering optical thickness 
of anv one layer docs not exceed 0.02. The discussion given above 
applies to a homogeneous atmospheric model, or to a model in which 
the nonhomogeneity varies slowly with height. If no' numerical 
results for a basic nonhomogeneous model divided into 160/n^ and 
160/^ number of layers may be found to differ significantly from each 
other as the division procedure used does not provide exactly identical 
models in both cases. 

( s r) 

The output for a non-absorbing, Rayleigh atmosphere (t, ' =0.1 

b 

and 1.0, 0=0 and 60 , R= 0.0, 0. 06941 and 0. 25) was checked 

o 

for several values of 6 and <p with that obtained using Chandrasekhar's 

X- and Y- function method by the author, and also by the Technical 

Officer of this contract, Dr. R.S. Fraser. A value of 2° was used for 

the parameter 66 in all cases. Values of the intensity of the radiation 

emerging at the top and at the bottom of the atmosphere were found to 

agree, on the average, within +3 units in the fourth place after the 

{s r) 

decimal point {~0.2% accuracy) when the models with t ' =0.1 

and 1.0 were divided into 16 and 40 layers, respectively. Values 
of the degree of polarization were found to agree within + 1 units in 
the third place after the decimal points, and those of the direction of 
polarization were found to agree within + 0.1° , 

Computations were also performed for a nonhomogeneous, 

Rayleigh case ( X = 0. 3 125fi , T^ S,r ^ = 1. 025, 6 q = 60° , total ozone = 0. 341 
cm -atm) for which numerical results of good quality are available in 
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published form (Dave and Furukawa, 19C6). For this run, the basic 
atmospheric model was divided into 80 layers, and a value of 2° was 
taken for the input parameter L h , Values of intensities obtained 
using the fourth vector program were found to agree, on the average, 
within +2% with those of Dave and Furukawa? and the values of degree 
of polarization within + 4 units in the third place after the decimal 
point. The results of this comparison are considered as satisfactory 
especially when one notes that Dave and Furukawa divided their 
model into 366 layers. However, they used a value of 1,000 mb for 
their surface pressure (versus 1,013 mb in our case), and they approxi- 
mated the contribution from higher- than-seventh-order scattered radia- 
tion by a geometric series. 
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VI. COMPUTER REQUIREMENTS 


In Table I we have listed some computing facilities required 
during test executions of the various programs. It may be added 
that a 9-track, 1600 BPI tape was used as an input or output tape 
during execution of each of the eight programs. The computer facility 
used for the test run was an IBM System/360 Model 91 computer 
located at the NASA Goddard Space Plight Center. 

During the execution of the first program, values of the Legendre 
coefficients and efficiency factors needed to fully describe the 
scattering properties of a single spherical particle were computed 
for the discrete values of the size parameter x given by x = 0. 1(0.2)169.9. 
A value of 1.5-0.031 was used for tht> .efractive index of the material 
of the particle. The CPU time needed for all computations for a given 
value of x increases with x; e.g. , CPU time for the vector case and 
for x=l.l, 10.1, 100.1, and 169.9 was 0.02, 0.09, 13.98, and 
57,73 seconds, respectively. For the scalar case, the corresponding 
timings were smaller by about 20%. 

During the execution of the second program, values of the 

normalized Legendre coefficients and absorption as well as scattering 

volume coefficients for a unit volume containing an arbitrary size 

distribution (Haze C; r . = 0.03u , r = 10. Ou, change of slope 

min max 

at 0 . lji , m = 1 . 5 - 0 . 03i) of spherical aerosols , were computed for 
monochromatic radiation of 0.55^ wavelength. The value of the 
parameter x corresponding to the largest particle is 114, and 
accordingly we required 267 terms in the Legendre series. 
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Table I 


Compiler, storage in kilobytes, use (or non-use) of 2314 direct 
access facility, and Central Processing Unit (CPU) time in minutes 
during test runs of various programs. 


Program 

Compiler 

Storage in 
kilobytes 

2314 direct 
access facility 

CPU time 
in min. 

Scalar I 

G 

200 

NO 

184 

Vector I 

G 

200 

NO 

224 

Scalar II 

G 

200 

NO 

0.4 

Vector II 

G 

200 

NO 

0.5 

Scalar III 

H,OPT=2 

220 

YES 

14 

Vector III 

H,OPT=2 

748 

NO 

20 

Scalar IV 

H,OPT=2 

SOO 

YES 

5 

Vector IV 

H,OPT=2 

750 

YES 

31 
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It may be pointed out that the first two programs require less 
than 20OK bytes of storage space during execution. A region of 200K 
bytes was asked for because the system used assigned a minimum 
of 200K bytes all the time. 

During the execution of the third program, values of the Fourier 
coefficients were computed for P’ = 0°(2 O )180° and 9=0°(2°)90° 
for the case for which punched card output was obtained after execution 
of the second program. Because of the use of a less strict criterion, 
this Fourier series output could be terminated after 225 terms instead 
of 267. 

For the test runs for the fourth program, values used for several 

input parameters of interest are as follows: 6 = 40° , 59=2° , 

(s r) ^ (s ml 

Number of layers = 40, X = 0.55m , = °' 098 ' T b = °* 242 ' 

(a , m) _ Q0Q ( no ozone, and a typical atmospheric aerosol number 
b 

density versus height curve with Haze C for the size distribution 

function. Because of the use of the modified Fourier series method 

[see discussion following Eq. (29)3, the number of terms of the series 

to be evaluated depends upon the nature of the aerosols and value of 

9 . For this particular case, we required 147 terms instead of 225. 

o 

The CPU time and the number of iterations required for evaluation of 
some selected terms of the Fourier series are given in Table II for 
the vector and scalar cases. 

Finally, it should be pointed out that this CPU time for the 
fourth program depends upon the number of layers into which the basic 
atmospheric model is divided. The above timings will be approximately 
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Table II 


CPU time and number of Iterations required for evaluation of 
some selected terms of the Fourier series during multiple-scattering 
calculations in a plane-parallel atmosphere (for details of the model, 
please see the text). 


Term of the 

CPU time in secs. 

Number of iterations 

series 

Scalar 

Vector 

Scalar 

Vector 

Lambert ground 

13.5 

108.9 

5 

5 

n = 1 

15.4 

123.3 

6 

6 

n - 10 

7.3 

53.2 

4 

4 

n = 50 

1.4 

6.8 

3 

3 

n = 100 

0.8 

2.3 

2 

2 

n = 147 

0.5 

1.6 

2 

2 


twice as much when calculations are repeated for an 80-layer model. 
On the other hand, these timing:: also depend very strongly on the 
value of the input parameter A 6 . The timings given above were 
found to decrease by a factor of about 8 and 20 when calculations for 
the above case were repeated for A 0 « 6° and 10°, respectiveiy. 
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